Why choose just one model for prediction?
fm1 <- brm(Energy ~ 1, data = M,
prior = prior(normal(0, 3), class = Intercept), iter = 2e4,
refresh = 0, seed = 3476283)
fm2 <- brm(Energy ~ Caste, data = M,
prior = prior(normal(0, 3), class = b), iter = 2e4,
refresh = 0, seed = 312379)
fm3 <- brm(Energy ~ Mass, data = M,
prior = prior(normal(0, 3), class = b), iter = 2e4,
refresh = 0, seed = 12365864)
fm4 <- brm(Energy ~ Mass + Caste, data = M,
prior = prior(normal(0, 3), class = b), iter = 2e4,
refresh = 0, seed = 8873542)
fm5 <- brm(Energy ~ Mass * Caste, data = M,
prior = prior(normal(0, 3), class = b), iter = 2e4,
refresh = 0, seed = 612356)mw1 <- model_weights(fm1, fm2, fm3, fm4, fm5, weights = "waic")
mw2 <- model_weights(fm1, fm2, fm3, fm4, fm5, weights = "loo")
mw3 <- model_weights(fm1, fm2, fm3, fm4, fm5, weights = "stacking")
mw4 <- model_weights(fm1, fm2, fm3, fm4, fm5, weights = "pseudobma",
BB = FALSE)
mw5 <- model_weights(fm1, fm2, fm3, fm4, fm5, weights = "pseudobma")
tibble(Model = paste0("fm", 1:5),
WAIC = mw1,
LOO = mw2,
Stacking = mw3,
PseudoBMA = mw4,
`PseudoBMA + BB` = mw5) |>
knitr::kable(digits = 2)| Model | WAIC | LOO | Stacking | PseudoBMA | PseudoBMA + BB |
|---|---|---|---|---|---|
| fm1 | 0.00 | 0.00 | 0.06 | 0.00 | 0.02 |
| fm2 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 |
| fm3 | 0.04 | 0.04 | 0.04 | 0.04 | 0.14 |
| fm4 | 0.50 | 0.50 | 0.56 | 0.50 | 0.42 |
| fm5 | 0.46 | 0.45 | 0.35 | 0.45 | 0.41 |
posterior_epred() (Lecture 3.4)posterior_predict()pp_fm4 <- crossing(Mass = seq(3.8, 5.3, length.out = 200),
Caste = levels(M$Caste))
pp <- posterior_predict(fm4, newdata = pp_fm4) |>
as.data.frame()
str(pp)'data.frame': 40000 obs. of 400 variables:
$ V1 : num 2.8 3.35 3.56 3.16 3.32 ...
$ V2 : num 3.36 3.76 3.5 4.27 3.13 ...
$ V3 : num 3.27 3.07 3.13 3.23 3.27 ...
$ V4 : num 2.93 3.85 3.34 3.57 4.18 ...
$ V5 : num 3.59 3.7 3.99 2.66 2.79 ...
$ V6 : num 3.69 3.59 3.06 4.02 3.44 ...
$ V7 : num 3.02 3.48 3.14 3.09 3.35 ...
$ V8 : num 3.98 3.85 3.92 4.25 3.38 ...
$ V9 : num 3.54 3.02 3.6 3.47 2.72 ...
$ V10 : num 4.02 3.74 2.85 4.05 3.07 ...
$ V11 : num 2.41 3.56 2.92 3.75 3.32 ...
$ V12 : num 3.68 3.72 3.74 3.97 3.36 ...
$ V13 : num 2.89 3.46 3.23 3.45 3.19 ...
$ V14 : num 4.29 3.79 3.72 3.9 3.72 ...
$ V15 : num 3.19 3.38 3.55 3.45 2.79 ...
$ V16 : num 3.78 4.15 3.12 3.75 3.39 ...
$ V17 : num 2.31 3.14 3.96 3.54 3.3 ...
$ V18 : num 3.75 3.92 4.08 4.46 3.88 ...
$ V19 : num 2.7 3.06 3.51 2.95 2.35 ...
$ V20 : num 3.43 4.02 3.89 3.86 3.73 ...
$ V21 : num 3.41 3.53 3.72 3.37 3.03 ...
$ V22 : num 3.69 3.67 3.52 3.63 4.12 ...
$ V23 : num 3.3 3.14 3.39 3.12 2.94 ...
$ V24 : num 3.79 3.84 3.86 4.13 3.34 ...
$ V25 : num 3.49 3.64 3.17 2.5 3.32 ...
$ V26 : num 3.11 4.01 3.67 3.99 3.65 ...
$ V27 : num 3.23 3.91 3.97 3.31 3.33 ...
$ V28 : num 3.01 3.48 3.67 3.89 3.71 ...
$ V29 : num 4.25 3.59 3.73 2.64 3.09 ...
$ V30 : num 3.89 3.29 3.46 3.49 3.34 ...
$ V31 : num 3.62 3.62 3.94 3.3 3.17 ...
$ V32 : num 3.39 3.6 3.83 4.01 3.94 ...
$ V33 : num 3.58 2.97 3.83 3.17 3.19 ...
$ V34 : num 3.31 3.8 3.57 4.32 2.99 ...
$ V35 : num 3.51 3.58 3.28 2.84 3.53 ...
$ V36 : num 3.95 3.34 3.5 3.86 3.94 ...
$ V37 : num 3.93 3.18 3.53 2.93 3.85 ...
$ V38 : num 3.35 3.61 3.8 3.74 3.39 ...
$ V39 : num 3.57 3.71 3.72 3.33 3.39 ...
$ V40 : num 3.94 4 3.89 4.49 3.82 ...
$ V41 : num 2.67 3.21 4.08 3.79 3.35 ...
$ V42 : num 3.55 3.19 3.29 3.62 3.5 ...
$ V43 : num 3.17 3.47 3.42 3.27 3.24 ...
$ V44 : num 3.83 3.97 3.16 4.1 4.09 ...
$ V45 : num 3.8 3.47 3.71 3.21 2.94 ...
$ V46 : num 3.48 3.98 4.25 3.62 3.97 ...
$ V47 : num 3.5 3.45 3.67 3.29 3.63 ...
$ V48 : num 3.68 4.03 3.85 3.8 3.97 ...
$ V49 : num 3.16 3.63 3.49 3.65 2.85 ...
$ V50 : num 3.3 4.61 3.85 3.37 3.85 ...
$ V51 : num 3.67 3.65 4.19 3.07 3.82 ...
$ V52 : num 3.9 3.75 3.53 3.83 3.96 ...
$ V53 : num 3.46 3.93 3.68 3.55 3.23 ...
$ V54 : num 4.22 4.01 2.88 4.04 4.03 ...
$ V55 : num 2.77 3.65 3.32 3.85 3.77 ...
$ V56 : num 3.55 4 3.26 3.2 3.52 ...
$ V57 : num 2.85 3.3 4.11 3.2 2.62 ...
$ V58 : num 3.04 3.84 3.68 4.01 3.94 ...
$ V59 : num 3 4.09 3.29 3.46 3.11 ...
$ V60 : num 3.73 3.88 4.2 3.82 4.03 ...
$ V61 : num 2.88 3.56 3.55 3.48 3.22 ...
$ V62 : num 4.26 4.02 3.46 3.45 4.09 ...
$ V63 : num 2.51 3.59 4.42 3.08 3.56 ...
$ V64 : num 3.74 3.88 3.96 4.1 3.38 ...
$ V65 : num 3.98 3.26 3.67 3.19 3.33 ...
$ V66 : num 3.52 3.8 3.86 3.94 3.92 ...
$ V67 : num 3.9 3.7 3.45 3.33 2.99 ...
$ V68 : num 4.4 3.84 3.97 3.83 4.13 ...
$ V69 : num 2.98 4.06 4.16 4.07 3.61 ...
$ V70 : num 4.34 4.4 3.49 4.07 3.35 ...
$ V71 : num 3.18 3.31 3.33 3.75 3.68 ...
$ V72 : num 3.09 3.46 3.37 3.98 3.39 ...
$ V73 : num 3.37 3.66 4.56 3.61 2.8 ...
$ V74 : num 4.06 4.1 3.43 4.38 3.5 ...
$ V75 : num 3.08 3.47 3.59 3.43 3.45 ...
$ V76 : num 4.6 4.14 3.65 4.14 3.85 ...
$ V77 : num 3.79 3.99 3.45 3.79 2.75 ...
$ V78 : num 3.5 3.44 4.5 4.07 4.07 ...
$ V79 : num 3.6 3.52 3.2 3.26 3.43 ...
$ V80 : num 3.25 4.45 3.95 3.74 3.7 ...
$ V81 : num 2.69 3.91 3.4 3.85 2.42 ...
$ V82 : num 3.99 3.96 3.66 3.74 3.51 ...
$ V83 : num 3.69 3.3 3.73 3.36 3.42 ...
$ V84 : num 4.18 4.09 3.8 3.9 3.5 ...
$ V85 : num 3.66 3.98 4.07 3.39 3.87 ...
$ V86 : num 3.69 3.84 3.84 4.24 3.93 ...
$ V87 : num 3.26 4.12 3.74 2.99 3.27 ...
$ V88 : num 3.45 3.83 4.43 4.19 3.85 ...
$ V89 : num 3.8 4.2 3.68 3.79 2.26 ...
$ V90 : num 3.57 3.45 4.03 3.96 3.88 ...
$ V91 : num 3.1 3.76 4.03 3.73 2.86 ...
$ V92 : num 4.17 4.2 3.94 4.22 3.7 ...
$ V93 : num 3.72 3.61 4.72 3.58 3.36 ...
$ V94 : num 3.56 4.35 3.89 4.18 3.64 ...
$ V95 : num 3.72 3.64 4.23 3.46 3.22 ...
$ V96 : num 4.86 3.86 3.53 4.28 3.65 ...
$ V97 : num 3.47 3.67 4.22 3.89 3.37 ...
$ V98 : num 3.75 4.22 3.5 4.09 3.85 ...
$ V99 : num 3.47 3.98 3.71 3.5 3.67 ...
[list output truncated]
pp_fm4 <- pp_fm4 |>
mutate(Q50 = apply(pp, MARGIN = 2, FUN = median),
Q5.5 = apply(pp, MARGIN = 2, FUN = quantile, prob = 0.055),
Q94.5 = apply(pp, MARGIN = 2, FUN = quantile, prob = 0.945))
head(pp_fm4)# A tibble: 6 × 5
Mass Caste Q50 Q5.5 Q94.5
<dbl> <chr> <dbl> <dbl> <dbl>
1 3.8 NonWorker 3.30 2.69 3.92
2 3.8 Worker 3.69 3.16 4.23
3 3.81 NonWorker 3.31 2.69 3.92
4 3.81 Worker 3.70 3.16 4.23
5 3.82 NonWorker 3.31 2.70 3.93
6 3.82 Worker 3.70 3.17 4.24
Combine all models by their relative support
pp_mod_avg <- crossing(Mass = seq(3.8, 5.3, length.out = 200),
Caste = levels(M$Caste))
pp <- pp_average(fm1, fm2, fm3, fm4, fm5,
newdata = pp_mod_avg,
weights = "stacking",
method = "posterior_predict",
probs = c(0.055, 0.5, 0.945))
attr(pp, "weights") |> round(3) fm1 fm2 fm3 fm4 fm5
0.057 0.000 0.038 0.556 0.350
fm1 fm2 fm3 fm4 fm5
2267 0 1521 22231 13981
# A tibble: 6 × 7
Mass Caste Estimate Est.Error Q5.5 Q50 Q94.5
<dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 3.8 NonWorker 3.45 0.473 2.73 3.43 4.25
2 3.8 Worker 3.70 0.364 3.14 3.70 4.30
3 3.81 NonWorker 3.45 0.467 2.74 3.43 4.24
4 3.81 Worker 3.71 0.363 3.15 3.70 4.30
5 3.82 NonWorker 3.46 0.470 2.75 3.44 4.26
6 3.82 Worker 3.72 0.362 3.16 3.71 4.31
Also Gelman et al. (2020) and Gabry et al. (Gabry et al. 2019)